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Abstract 

The Vlasov equation is a kinetic model describing the evolution of a plasma which is a 
globally neutral gas of charged particles. It is self-consistently coupled with Poisson's equation, 
which rules the evolution of the electric field. In this paper, we introduce a new class of forward 
Semi-Lagrangian schemes for the Vlasov-Poisson system based on a Cauchy Kovalevsky (CK) 
procedure for the numerical solution of the characteristic curves. Exact conservation properties 
of the first moments of the distribution function for the schemes are derived and a convergence 
study is performed that applies as well for the CK scheme as for a more classical Verlet scheme. 
A L 1 convergence of the schemes will be proved. Error estimates (in 0(At 2 + h 2 + ^) for 
Verlet) are obtained, where At and h = max(Ax, Av) are the discretisation parameters. 
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1 Introduction 

The Vlasov equation describes the dynamics of charged particles in a plasma or in a propagating 
beam. The unknown f(t, x, v) which depends on the time t, the space x and the velocity v represents 
the distribution function of the studied particles. The coupling with the self-consistent electric fields 
is taken into account through the Poisson equation. 

The numerical solution of such systems is most of the time performed using Particle In Cell (PIC) 
methods, in which the plasma is approximated by macro-particles (see [1]). They are advanced in 
time with the electromagnetic fields which are computed on a grid. However, despite their capability 
to treat complex problems, PIC methods are inherently noisy, which becomes problematic when 
low density or highly turbulent regions are studied. Hence, numerical methods which discretize 
the Vlasov equation on a grid of the phase space can offer a good alternative to PIC methods 
(see [U dU [T21 [T71 [5] ) . The so-called Eulerian methods can deal with strongly nonlinear processes 
without additional complexity, and are well suited for parallel computation (see [14J). Moreover, 
semi-Lagrangian methods which have first been introduced in meteorology (see |16} [I~8| I19j). try 
to take advantage of both Lagrangian and Eulerian approaches. Indeed, they allow a relatively 
accurate description of the phase space using a fixed mesh and avoid traditional step size restriction 
using the invariance of the distribution function along the trajectories. 

Traditional semi-Lagrangian schemes follow the characteristics backward in time. In [7], follow- 
ing the idea of Reich [15] , we introduced a forward Semi-Lagrangian scheme for the Vlasov- Poisson 
system based on a forward numerical solution of the characteristics using a classical Verlet or 
Runge-Kutta (order 2 and 4) scheme. The Verlet scheme can only be applied for specific differen- 
tial equations, as for example the characteristics of the Vlasov-Poisson system, but not for more 
general cases, as the characteristics of the guiding centre or the gyrokinetic approximation of the 
Vlasov equation. Therefore an alternative to Verlet is necessary. On the other hand, Runge-Kutta 
schemes, which can be used in the general case, are very costly in our context, especially when 
going to higher order, as they require a deposition of the charge and the solution of the Poisson 
equation at intermediate time steps. We propose here, a new scheme for the characteristics based 
on a Cauchy-Kowalevsky (CK) procedure, that can be performed up to an arbitrary order. Second 
and third order are developed in the present paper. We shall also discuss the conservation of the 
first moment for both Verlet and CK algorithms. 

A proof of the convergence of PIC method for the Vlasov-Poisson system was performed by 
Cottet and Raviart [9]. Proofs of convergence and stability of the classical Semi-Lagrangian method 
applied to the same model were obtained by Besse and Mehrenberger [I]. These estimates are made 
in L 2 norm, since L°° seems out of reach as they explain. They manage to do it because they deal 
with split methods, and thus only consider constant coefficient transport at each split step. In order 
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to prove convergence in more general cases, the L norm seems appropriate, as it enables to use 
the partition of unity property of the splines. Moreover, Despres [TO] explains possible advantages 
of studying L 1 convergence instead of more common L 2 . 

We propose here a proof of L 1 convergence of the forward semi-Lagrangian scheme with both 
Verlet and CK solution of the characteristics in the particular case of linear spline interpolation. 
We also obtain second order error estimates in time and space. 

This paper is organized as follows. In the first part, the continuous problem is presented. In 
the second part, the discrete problem and the numerical scheme to solve it are explained. We also 
prove the exact conservation of the first moment with respect to v at the discrete level for both 
CK and Verlet schemes. Then the convergence of our numerical schemes is proved and finally the 
schemes are validated and compared on a couple of classical test problems. 



2 The continuous problem 

2.1 The Vlasov-Poisson model 

Let us consider f(t, x,v) > the distribution function of positively charged particles in phase-space, 
and E(t, x) the self consistent electric field. The dimensionless Vlasov Poisson system reads 

df 

_± + vdxf + E (t,x)d v f = 0, (2.1) 



8 x E(t,x) = p(t,x)= [ f(t,x,v)dv-l, (2.2) 



where x and v are the phase space independent variables. A periodic plasma of period L is con- 
sidered. So x G [0, L], v € R, t > 0. The functions / and E are submitted to the following 
conditions 

f(t,0,v) = f(t,L,v),Vv€R,t>0, (2.3) 

E(t, 0) = E(t, L) o ~ [ f f(t,x, v)dvdx = l,Vi > 0, (2.4) 
L Jo Jr 

which translates the global neutrality of the plasma. In order to get a well-posed problem, a zero- 
mean electrostatic condition has to be added, which corresponds to a periodic electric potential: 



f E(t,x)dx = 0, W>0, 



(2.5) 



and an initial condition 

f{0,x,v) = f (x,v), Vxe[0,L],veR. (2.6) 

Assuming that the electric field is smooth enough, equations (|2.ip . (I2.3P and (12. 6j) can be solved in 
the classical sense as follows. 

The first order differential system 

dX 

— (t;(x,v),s) = V(t;(x,v),s), 
dV 

— (t;(x,v),s) = E(t,X(t;(x,v),s)), (2.7) 

where (X(t;(x,v),s),V(t;(x,v),s)) are the characteristic curves, solutions of (|2.7|) at time t with the 
initial condition 

X(s; (x, v), s) = x, V(s; (x, v),s) = v. (2-8) 
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For the existence, the uniqueness and the regularity of the solutions of this differential system, the 
reader is referred to R| • The solution of problem (|2.ip . (|2.6p is then given by 

f(t,x,v) = f (X(0;(x,v),t),V(0;(x,v),t)), Vx G [0, L},v G R, t > 0. (2.9) 

Since 

d(X,V) _ 
d(x, v) 

the conservation of particles is ensured for all times: 

— I J f(t,x,v)dvdx = — I J fo(x,v) dv dx = 1. 
L Jo Jr L Jo Jul 

According to previous considerations, an equivalent form of the Vlasov-Poisson periodic problem 
is to find (/, E), smooth enough, periodic with respect to x, with period L, and solving the equations 
(2.2), (|2.7p . (|2.8p and (|2.9p . Introducing the electrostatic potential ip = (p(t,x) such that E(t,x) = 
— d x (p(t,x), and setting G = G(x,y) the fundamental solution of the Laplacian operator in one 
dimension. That is —d x G(x,y) = 6q(x — y) with periodic boundary conditions. It comes 

E(t,x)= [ K(x,y)(f f(t,y,v)dv-l)dy, 
JO JR 

where 

K(x,y) = -d x G(x,y) = (| - 1), < x < y, 

= y, y < x < L. 

2.2 Existence, uniqueness and regularity of the solution of the continuous prob- 
lem 

Theorem 1 Assuming that fo G Wc^per^ (Re x R>„ 

) ^ per^ (Ra; x IR^ being the Sobolev space of 
functions with first derivatives in L°° , compactly supported in v and periodic in x), positive, periodic 
with respect to the variable x with period L, and Q(0) < R, with R > defined as follows 

Q(t) := l + sup\v\ : x G [0,L],t G [0,t]|/(r,z,u) ^ 0, 

and 

- / fo{x,v) dvdx = 1, 

then the periodic Vlasov-Poisson system has a unique classical solution (f,E), periodic in x, with 
period L, for all t in [0, T], such that 

E g W 1,oc (Q, T; Wp£?(M)), 
and there exists a constant C = C(R, /o) dependent of R and /o such that 

Q(T) < CT. 

Moreover if we assume that f G W™ P ^ X (R X x R v ), then (f,E) G W m >°°(0,T;W™ p ^ x (R x x R v )) x 
W 1 ' 00 ^, T; Wp£° x (R)), for all finite time T. 

For the proof, the reader is referred to the references [13] . 
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3 The discrete problem 



3.1 Definitions and notations 

Let Q = [0,L[x[— R, R], with R > Q(T), and a cartesian mesh of the phase-space VI. Mh is 
given by a first increasing sequence (aJi)ie[o..jvJ of the interval [0,L] and a second one (vj)je[o..N v ] 
of the interval [-R,R]. Let Ax^ = Xj+i — Xi the physical space cell width and Avj = Vj+i — Vj 
the velocity space cell width. In order to simplify the study, a regular mesh will be used, i.e 
Axi = Ax = N L +1 , and Avj = Av = j^, where N x , N v belong to N. Then h is defined being 
max(A:c, Av). 

For each function g defined on all the points (xi,Vj) 6 we will set gij := g(xi,Vj), and the 
sequence is completed on Z x Z by periodicity in x and by in v. The sequence (xi,Vj) will also 
be defined on the whole set Z x Z, by Xi := iAx, and Vj := —R + jAv. The set of all L-periodic 
functions in x and compactly supported in v will be denoted P(Vl). 

Now let / = (/i,j)(ij)eZxZ be a continuous grid- function, periodic in the x direction and 
compactly supported in the v direction, with a support included in [—R,R]. If / is a function 
defined on the points M^, a discrete grid-function / can be defined by /y := f(xi,Vj) for all 
(i, j) G [0, iVa.] x [0, N v ]. In order to lighten notations, / will be kept instead of /. Let L%(Q) (resp. 
L\(Q)), the set of grid-functions whose 11-11^2^^ (resp ||.||£i(Q)) is bounded 

N x N v 

\\f\\ Lm = (AxAvJ2Y,\hi\ 2 ) h > 

i=0 j=0 

N x N v 

\\f\\ Li{Q) = AxAvJ2Y,\hi\- 

i=0 j=0 

As was precised in the introduction, a L 2 convergence analysis for backward Semi-Lagrangian 
scheme, in the case of a Strang split time advance, was performed in pp. In paper [7], it is explained 
that for split methods, where the split steps consist of constant coefficient transport, forward and 
backward methods are exactly the same. So all the L 2 results exposed in [Tj are also valid for our 
method when time splitting is used. In this paper, we shall consider the convergence of a non split 
method, and the L 1 norm seems more appropriate for this kind of study. 

Remark 3.1 If f € L|(J2), then f G L\(Vl) thanks to the Cauchy Schwarz inequality, as Vl is 
bounded. 

In the sequel, a final time T is fixed, as well as a uniform time discretization {t n ) n <N T of the 
interval [0, T], with time step At = t n+1 — t n . At each point (xi,Vj) G M^, an approximation 
fh(t n ,Xi,Vj) of the exact distribution function f(t n ,Xi,Vj) at time t n = nAt is defined. The 
approximation function fh(t n ) is then given at each point of ~R X x M v thanks to an interpolation 
operator Rh defined on a uniform grid: 

R h : L^Vl) n P(fi) — ► L x (fi) n P(Q), 

f^R h f= /<j*<j. 

(i,j)ezxz 

where VPj will be linear spline functions for our study. In numerical results, since linear interpola- 
tion is quite diffusive, cubic splines will be used. In order to get a convergent scheme, the operator 
Rh must satisfy some approximation properties which will be detailed later. 
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3.2 The numerical scheme 

The electric field operator for the real- valued function g £ -^ 1 ([0, L) x R) is defined this way: 

E[g](x) = [ L K(x, y){ [ g(y, v)dv - 1). (3.10) 

The approximate function is solution on the grid of the following Vlasov equation: 

^(i, x, v) + v -^(ti x > v ) + E h(t, x)^-(t, x, v) = 0. 
This function follows approximate characteristics, solutions of 

-^(*; ( x > s ) = v k{t\ ( x , v)i s )> 

-J^(t; (x, v), a) = £ h (i, X(i; (x, t;), a)), (3.11) 

where is defined exactly from fh using (|3.10p : Eh = E[fh](x). So we get: 

Vi E [i n , t n+1 [: / fc (t, x, «) = ^ <^(x - (x fc , t;,), t n ))5 h (x; - V h (t, (x k , Vl ),t n )), 

k,i 

so that is given on the mesh at time t n by: 

fh(t n ,Xi,Vj) = ^2u]% tl S h (xi - x k )S h (vj - vi) Vn. 
k,l 

These are the interpolation conditions enabling to define fh everywhere. The computation of 
( u kl)kl from the grid values amounts to solving a linear system, which is trivial in the case of 
linear splines, where u; k l l = fh(t n ,x k ,vi). 

Let us recall that the linear B-spline S is defined as follows 

s(x) = l (1 " N) tfo^N^ 1 . 

1 otherwise. 
The spline Sh actually used, which shall be defined with the size of the mesh will be 

Sh(x) = 

with h = Ax for splines in the x variable and h = Av for splines in the v variable. From now on, 
Sh will be denoted by S for the sake of simplicity. 

The distribution function is updated this way: The ending point of the characteristic starting 
from (xi,Vj) is computed: (Xh(t n+ ; (xi,Vj),t n ), V/ l (i n+1 ; (xi,Vj),t n )),\/(i, j). Then, since fh is 
constant along the approximate characteristics, the value is deposited on the nearest grid points, 
the number of which depending on the degree of the splines used for the interpolation. This amounts 
to computing fh at time i n +i at the grid points using the following formula. 

fh{t n+ \x hVj ) = Y, fh(t n , x k , Vl )S( Xi - X h (t n+1 ; (x k , Vl ),t n ))S{ Vj - V h (t n+1 ; (x k , Vl ),t n )) V(i, j), 
k,l 

Note that 

(X /t (t"+ 1 ;(x i ,^-),i n ),^(t" +1 ;(^,^),i n )),V( i ,i), 
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are computed by a numerical solution of the differential system (13. IIP , Since this requires an 
explicit solution of that system, any standard ODE solver such as Verlet, Runge-Kutta or others 
can be used. Our analysis will be based on the Verlet algorithm, which is second order accurate, 
and on a Cauchy Kovalevsky procedure, which can be of any order, as an alternative to more costly 
Runge-Kutta solvers. But we will consider only the second and third order, because higher ones 
would not increase accuracy in our case, as we will explain. 

3.2.1 Verlet algorithm 

Starting at time t n from the grid point (2^,1^) 

• Step 1: Vfc,Z, xj!j~ 2 - Xk = ^fVi, 

• Step 2: compute the electric field at time t n+ z 

— deposition of the particles on the spatial grid xi for the density p^: ph{xi,t n+ i) = 

J2k i w fc i^\ x i ~ x k i 2 )' n ^ e ™ a PIC method. 

— solve the Poisson equation on the grid xf E(xi,t nJr i). 

• Step 3: VM, vffl-vi = At E(x^ V n+ 5), 

. Step 4: Vk,l, xlf - x n k f = f vft\ 

This is the way the algorithm is implemented. In our convergence study, the slight difference is 
that an exact solution of Poisson's equation, based on the Green formula (|3.10p . is used. 



3.2.2 Cauchy Kovalevsky procedure 

The idea is to get high order approximations of the characteristics using Taylor expansions in time. 
And then, using the charge conservation equation, and higher velocity moments of the Vlasov 
equation, to replace time derivatives with terms containing only spatial derivatives and moments 
at time t n which can be easily computed. Up to third order, these Taylor expansions in time lead 
to 

At 2 At 3 d 

X n + 1 = X n + At yn + E n {x n )+ E{x{tht) 

2 b at 1 

At 2 d At 3 d 2 

yn+l = r + A ^„ (r) + __ £(I(t))t)|i=j „ + T# £(I(f),0 



t=t n 



In order to be able to compute all terms of these expansions we need the three first total time 
derivatives of E(X(t),t). 

jE(X{tlt) = ^(X(t),t) + ^(t)^(X(t),t) 

= -J(X(t),t) + J(t) + V(t)p(X(t),t), 

where p(x,t) = J f(x,v,t)dv — 1, J(x,t) = J f(x,v,t)vdv and J(t) = ^ Jq J(x, t) dx. Indeed, 
the Poisson's equation yields ^| = p and integrating the Vlasov equation with respect to velocity, 
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yields the charge conservation equation §? + = 0- Hence taking the derivative of the Poisson's 
equation with respect to time and using this equation we get 

9 , dE n n 

From which we obtain, as E(x,t) dx = 0, that 

dE 

The second order total derivative in time of E reads 

+ E(X(t),t)p(X(t),t) + V(t)(^(X(t),t) + V(t)^(X(t),t)). 

In order to use this expression, we need ^(t). 

The Cauchy-Kovalevsky procedure consists in getting rid of time derivatives, replacing them 
with space derivatives obtained from the equation, in our case, we use the velocity moments of the 
Vlasov equation. First for p, we use the charge conservation equation: 

^(X{t),t)=~(X(t),t). (3.12) 

In order to get the time derivative of the current J, we need to use the Vlasov equation, multiply 
it with v, and integrate it with respect to v, so that we get: 

dJ d r „ f df , 

— + —I 2 + E jfvdv = 0, 

at dx J R dv 

where I n (x,t) = J R f(x,v,t)v n dv so that, using that / is compactly supported and integrating by 
parts: 

— (X(t), t) = --Jl(X(t),t) + E(X(t),t)(l + p(X(t),t)). (3.13) 
Let us prove that ^ (t) = 



dJ, v 1 f L dJ, . , 

-d~t [t) = lJ m {t > x)dx 





1 I ^ dl 

= - (-^{t,x)+E(t,x)(l + p(t,x))dx 

= y([/ 2 (t,0) -I 2 (t,L)} + f E(t,x)dx+ [ E(t,x)p(t,x)dx 
L Jo Jo 

= j[E\t,L)-E 2 (t,0)] 

= 0. 

thanks to periodicity, in fact (|2.3|) . (|2.5|) . We will see later, that numerically this value is also zero. 
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We finally get the following third order Cauchy Kovalevsky (CK3) time algorithm, using (|3.12|) . 

At 2 At 3 

x n+l = X n + M yn + J_ E n, x n- ) + (V n p n (X n ) - J n (X n ) + J), 

2 6 
V n+i = v n + AtE n iX n ) + -^-(V n p n (X n )-J n (X n ) + J) 

A 4.3 ar a t a n 

+ (^(X n ,t n ) - E n {X n ) - 2V n ^-(X n ,t n ) + (V n ) 2 ^f-(X n ,t n )). 

6 ox ox ox 

Let us now introduce a notation, which will be useful later: 

A+2 A f 3 
X n+1 =X n + At yn + ^^(jn) + ^^(J^y^ 

2 6 

and 

At 2 At 3 

yn+l = yn + AtE n< X n^ + 0»(X n , V n ) + V n ), 

2 6 

where 0, 99 are naturally defined. 

Remark 3.2 Obviously, in order to get a second order algorithm (CK2), we just keep the terms 
until At 2 included. 

3.3 Exact conservation of number of particles and momentum 
3.3.1 B-spline interpolation 

First, let us recall some useful properties of B-splines interpolation. The linear space of B-splines 
of order m + 1 writes, denoting by s^ the m th derivative of s 

S m +i,Ax = i s ( x ) e C m -HR), S ( m+1 )(x) = 0,Vx G ( Xi ,x i+1 ),Vi G R}, 

if m + 1 is even, and 

Sm+l.A, = {s(x) G r'HK), S ( m+1 )(x) = 0,\fx G (x._i,x i+ i),Vf G R}, 

if m + 1 is odd. 

The space of B-spline functions in two dimensions is defined as the tensor product of ID spaces. 
Let us precise the interpolation operator: 



R h i:j (f)(x,v) = uJij(f)S(x - Xi)S(v - Vj), 
R h f(x,v) = Y,Rh l Jf)(x,v). 



h3 

Now come the properties: 



§m+l,h = Span(S m+ i(. - Xi)S m+ i(. - Vj);V(i,j) G Z), 
§ m +i,h C W k ' p 1 < p < 00 0<k<m, 

Stability I \RhfWvw <C\\f\\ LP(n) V/G i7(l])n?(l]), 1 < p < 00 (i), 
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• Consistency and accuracy. There exists C > j ||/— Rhf\\w k <p(n) — Ch m+1 k \f\w m + 1 -p(fi) Vf € 
W m+l,p(ty n p(ty i^p^oo o<fe<m (ii), 

• J2i s m(- - Xi) = 1 (iii), / S m (u)du = Ax (iv). 

• T,i v l s i( v l -v) = v. (v) 

For the last item, we will give the proof: Let us suppose v = v p + aAv, a £ [0, 1[ 

viSi(vi — v) = VpSi(aAv) + v p+ iSi(Av - aAv) 

i 

= (p + aAv) 
= v 

Let us also precise particle and momentum conservation. The proof for the mass is independent 
from the spline degree, and the one for the first moment will only be shown for linear splines, even 
though it has been checked for the first three splines. 

3.3.2 Particle conservation 

The discrete algorithm preserves the total number of particles, as the following computation shows: 

,n+l 



m 



= J f h (t n+1 ,x,v)dxdv, 
= ^I*/ 1 / S(. x ~ x i)S( v ~ v j) dx dv, 
= AxAv^fh^i^vj), 



= AxAv £ "k,iS(xi - X(t n+1 ; (x k , v t ), t n )S{v - V(t n+1 ; (x k , Vl ),t n ), 

i,j k,l 

= AxAvY J ^k,i = ^AvY J f n {x i ,v j ) 



m 

k,l i,j 



thanks to partition of unity property (iii). 

Let us precise the way the interpolation operator acts, in fact: 



f h (t n+1 ,x,v) = R^<j S ( x ~ X h(t n+1 -,(x llVj ),t n )S(v -V h (t n+1 -,(x t , Vj ),t n ), 

= ^u^Six-x^v-vj), ( 3 - 14 ) 



h3 

by definition of ufj. This implies a kind of continuity of fh at time t n on the grid points. 
3.3.3 Momentum conservation 

Let us precise that in this paragraph (Xh(t; (xi,Vj),t n ),Vh(t; (xi,Vj),t n )) will be denoted (Xij(t), Vij(t)), 
and that Poisson will not be solved exactly. The aim here is to prove that Vn 



£ Vj f h (t n , x u Vj ) = v,f h (t n+ \x u Vj ). (3.15) 



1,3 h3 
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Let us distinguish two phases: the transport one and the deposition one. Let us start with the 
deposition phase, where we have to get: 

wg-^jSiCz* - xOSifa - vj) = Yl <j v i S i( x k ~ XiAt n+1 ))Si(vi ~ V i>j (t n+1 )) 

i,j,k,l i,j,k,l 

E ^j l viSi{x k - Xi)Si(vi - Vj) = ^ uj^ 1 ^ ViSi(vi - vj) 

i,j,k,l i,j I 

= E</S 

thanks to the property (v) of linear splines. Moreover 

*£v l S 1 (v l -V i>j (t n+1 )) = V i , j (t n+1 ), 
I 

thanks to the same property. So we finally get for the deposition phase: 

E< + S = ^<^(t n+1 ;(x,,^),t"). (3.16) 

Remark 3.3 This proof is given for linear splines, but was also checked for quadratic and cubic 
ones. The transport phase is independent from the spline degree. 

There remains to prove that 

»::,<•, = E<i y (* n+1 ; o*. ( 3 - 17 ) 

which corresponds to the transport phase. Note that this phase exists also in PIC methods, and 
the following proof of conservation of moments is adapted from [1] . 

Verlet We have with our Verlet algorithm: 

V(t n+1 ;(x l ,v j ),t n ) = v 3 +AtE n+1 2(X(t n+l 2,(x i ,v j ),t n )). (3.18) 

The electric field is only known on the mesh. In order to know it everywhere, we use a convolution 
between a spline function and the discrete E. 

E n+\ ( X ( t n+i f {Xij Vj)i r)) = E n +h (x k )S(x k - X(t n+ ^, ( Xi , Vj), t n )). 

k 

To get (I3.17D using (13.180 we just have to prove that 

]T ^E n+1 2 (x k )S(x k - X(t n+ K ( Xij Vj ),t n )) = 0, 

i,j,k 

and 

Y,^E n+ Hx k )S(x k -X(t n+ h } ( Xi , Vj ),t n )) = Ys^Hxk^^Sixk-Xit^^ixuVj),^)), 

i,j,k k i,j 

= ^2E n+ Hxk) P n+ hxk), 

= o, 
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for most of the centered algorithms used to solve Poisson numerically, like the following centered 
finite difference one on staggered mesh, with linear regularization: 

E n+ ^ (x k+ i) - E n+1 i (x k _i) = Axp n+ ? (x k ) Vfc, 

and 

E n+1 Hx) = Y,E n+ Hx l )S 1 (x - x t ). 

i 

Indeed, we get: 

^E n+ Hxk)p n+ Hx k ) = ^^^(^(^((S!^ -x i )-S 1 (x l ^i - Xi )), 

k k i 

i 

= li^^k^E^Hx^-^E^^E^Hx^)), 

i i 

= o, 

thanks to periodicity. 

To conclude, using (|3,16p and ()3.17p . we get (|3.15p . which is what was wanted. 

CK algorithm. We still have to prove that: 

That means for the third order scheme: 

Y,^,E n { Xl ) = (i) 

^2uijip n (xi,Vj) = (m) 

^2<^i,j4> n (xi,Vj) = (in) 

each number being linked with the order of the algorithm. 

First order. Using the same strategy (regularization of the electric field and centered algorithm): 

^2^E n (x t ) = J2^E n (x k )S(x t -x k ), 

i,j i,j,k,l 

= ^E n (x k )p n (x k ) = 0. (3.19) 
k,l 
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Second order. For the second order, we need J™,pf, J 

pi = Av^^k,iS{xi-x k )-l, 



k,l 



AxA 



J™ - J = Av Y ^k,iViS(xi - x k ) — Y ul t l v h 

k,l k,l 

Y j (>"i-r,) = Av uj^VjSixi - x h ) - Y ufyvj, (3.20) 

Y^UijJ n {xi) = Y Ui,j"k,lVlS(xi - X k ) - Yl W *>i>'> 

= AvY<^k,mS( Xl -x k )-Y^ (3-21) 

i,j,k,l k,l 

using mass conservation and Ylij^ij = L Q3.20n and (I3.21D are the same, just exchanging 
and (k, I). So 

Third order. Here we need I^ixi) 

h(t n ,Xi) = / v 2 o;g t S (x k - Xj)S(v - vjjdv, 
7ffi k,i 

= J2 wj? fifa - Xi ) f (v 2 S(v) + 2vviS{v) + v 2 S(v)) dv, 
k,i 

= Y^iSixk-Xijia + Avvf). 
k,l 

In <fr we still have Y2% j ^fjE n (xi) = 0. We also have three terms in -jj- which will be approached 
with a centered finite difference formula: 



'd x K ' ' 2Ax 

= 2A~x~ S W "X»( Q + ^X-SWi - x fc ) - S(**-i - x fc )), (3.22) 

i,j,k,l 

r) T 1 

2 E<^' 3-^) = Ai E - x fc ) - 5(x 4 „! - a*)), (3.23) 

Y^jv^it^xt) = ^Y <^k,iv](S(x l+1 - x k ) - S(x f _! - x fc )). (3.24) 
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Adding (pT22]h flH55j) and (1^241) and using (^T9]l we have: 

^2^i,j<t> n ( x i) = u iJ u k,l( a + v ? + 2w i^ + v j)( s ( x i+i ~ x k) ~ S(Xi-i - x k )), 

i,j i,j,k,l 

Yl w & w Jw( Q + u ? + 2u i^ + v j) s i x i+i ~ x k) = Yl u k<t u ij( a + v ? + 2w i v / + Vj)S(x k+1 - Xi) 

i,j,k,l ij,k,l 

= Y ^tj^kM + v i + 2v j v i + v j)S(x k - a?i_i), 

i,j,k,l 

just changing and (A;, /) and ^(xfc+i — Xj) = S^Xfc — Xj_i). So we get: 

Remark 3.4 We can see i/iai i/ie conservation of the first moment in v implies that numerically 
^ = 0, which means that ^ x ^ v UJ ?,j v j * s constant. 

4 Convergence analysis 

Theorem 2 Assume that /q G W^^ ra . (M x x positive, periodic with respect to the variable x, 
with period L, and compactly supported in velocity. 

Then the numerical solution of the Vlasov Poisson system (fh,Eh), computed by the numeri- 
cal scheme introduced in section 3.2 converges towards the solution (f,E) of the periodic Vlasov- 
Poisson system, and there exists a constant C = C(||/||vi/i,oo(o,T;VF 2 ' o °(f2))) independent of At and 
h such that for Verlet and CK2 algorithms: 

h 2 

11/ - fh\\l<x>(a,T;lMQ)) + \\E — Eh\\l°°(p,T;L<=°([0,L])) ^ C(At +h + -^). 
For CK3, we have: 

h 2 

11/ - fh\\l°°{0,T;lMp)) + \\ E ~ E h\\l°°(0,T;L°°([0,L})) < C(At +h + —). 



Remark 4-1 In order to get these estimates for CK, we will have to assume At < Ax. 
4.1 Decomposition of the error 

Let / be the exact solution of the Vlasov Poisson equation and fh the approximate solution pre- 
viously defined. In order to apply a discrete Gronwall inequality we express the I 1 error at time 

t n+l 

e n+ \i,j) = \f(t n+1 ,Xi,Vj) - f h (t n+1 ,Xi,Vj)\ \f(i,j), 

e n+1 = AxAv^V +1 (i,j). 

hi 
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Then f(t n+1 ,x k ,vi) — fh(t n+1 ,x k ,vi) can be decomposed as 



f{t n+ \x k , Vl )- f h {t n+ \x k , Vl ) = f(t n+ \x k ,v l )-R h f(t n+1 ,x k ,v l ) + 

R h f(t n+ \x k , Vl ) - R h f h (t n+1 ,x k , Vl ) + 
Rhfh(t n+ \x k , Vl ) - R h f h (t n+1 ,x k , Vl ), 



(4.25) 



where fh is the function fh at time t n but then follows the exact characteristics. Since f£ already 
belongs to the image of Rh, we have Rhfh(t n+1 -> x k-, v i) = fh{t n+l i x ki v l)- 

In order to estimate e n+1 , the three terms of the right hand side of the previous equation have 
to be dealt with. These estimations are developed in the following subsection. 

4.2 A priori estimates 

4.2.1 Stability for linear splines 

Let us translate the useful spline properties in this case, and give a few more results about the 
operator R^. 

Lemma 1 The Rh operator is consistent, that is, using property (i), for 1 < p < oo, and < k < 1 



This result is a classical property of B-splines. 

Lemma 2 With linear splines, if tOij(fo) > V(i,j) then\/n ujij(f n ) > \/(i,j). 

Proof: With linear interpolation, we get in fact Uij(f n ) = f n (xi,Vj), so if /o is positive, Uij(fo) 
is also, and since f n+1 (xi,Vj) is a sum of positive contributions coming from the f n (x k ,vi) which 
are positive by a recurrence hypothesis, it will also be positive, and so Uij(f n+1 ) is positive for all 
(i, j), and recurrently for all n. 

Lemma 3 Stability: Let f belong to C{Vi) n P(Q), then we have: 



3C>0 | \\f-R h f\\ wk , P{Q) <Ch 2 - k \f\ W 2, P{n) V/e^(fi)nP(fi). 



\\Rhf\\mn) = \\f\\ L l(ny 



Proof: 




— Il/llz,i(n)> 

using J S(x)dx = Ax, the positivity of / thanks to Lemma 2 and the positivity of /q. 
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4.2.2 Towards Gronwall 

Let us precise that in this subsection, some lemmas are valid for all the time algorithms we use, and 
when they are not, the lemmas will be proved in each case successively. For the Cauchy Kovalevsky 
procedure, the proofs will be done for (CK 3), since their adaptation to lower orders is trivial. We 
will now give estimates about the three right-hand side terms of the error e n+1 (|4,25p : 

Lemma 4 Let f belong to C(Q) n P(£l), then we have: 

\\f-R h f\\ LU n)<Ch 2 . (4.26) 

Proof: Thanks to Lemma 3 

H/-^/llLi(n) = WMf ~ Rhf)\\mn), 

< C\\f - R h f\\ LHQ ), 

< C / (ll/IU«(o,T;W 2 .° o (fi))^ 2 ! 

thanks to the property (ii) of spline interpolation and the fact that the domain is bounded. 
Lemma 5 Let f belong to C(fl) R P(fi), then we have: 

\\Rhf n+1 -R h fh n+1 \\ Lii u)<e n . (4.27) 

Proof: We compute 

\\R h f n+1 - Rj h n+ \ iiu) =AxAvJ2\(Rhf n+1 -Rjh n+1 )(x k , Vl )\, 



k,l 

= AxAv ^2 I y~](f n+1 ( x i> v i) ~ fh n+1 (xi,Vj))S(x k - Xi)S(vi - vj)\, 
k,l i,j 

i,j k,l 

< AxAvJ2Mf n )-"iAfh)l 

< e n , 

thanks once more to the partition of unity (hi) and f(xi,Vj,t n ) = Wjj(/ n ). 
Lemma 6 Let f belong to C(fl) D P(£l), then we have: 

\\R h fT +1 - R h f h n+1 \\ L{m < Cm^(\X(t n+1 ;(x l ,v 3 ),t n ) - X h (t n+1 ;(x t , Vj ),t n )\, 

hK ' 1,3 

+\V(t n+1 ;(x l ,v 3 ),t n )-V h (t n+1 ;(x u v J ),n\). (4.28) 

Proof: 



\Rj h n+1 - RhfJt +1 \\ L Un) = AxAvY,\(RJh n+1 - R h f h n+1 )(x k , Vl )\, 



AxAv I Yl "iAh) n {S{x k - X(t n+1 ; {x uV] ),t n ))S{vi - V(t n+1 ; (x u Vj ),t n )), 

k,l i,j 

S(x k - X h {t n+l - (a*, t n ))S(vj - V h {t n+l - ( Xi , Vj),t n )))\. 
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We can rewrite 



(S(x k - X(t n+1 ; (xuv&rySfa - V(t n+1 ; (x t ,v,),t n )) 

- S(x k - X h (t n+1 ; (x uVj )^))S(vi - V h (t n+1 ; {x uVj ),t n ))) 
= (S(x k - X(t n+1 ; ( Xi , Vj ),t n )) - S(x k - X h (t n+1 ; (x h Vj ), t n )))S( Vl - V(t n+1 ; (x hVj ),t n )) 
- (S( Vl - V(t n+1 ; (x uVj ),t n )) - S(v t - V h (t n+1 ; (*»,«,•), * n )))S(a* - X h (t n+1 ; (x uVj ),t n )). 

Then, we use the fact that Si is 1-Lipschitzian, compactly supported, and the property (i): 

\(S(x k - X(t n+1 ; (xi, Vj ),t n )) - S(x k - X h (t n+1 ; (x h Vj ), t n )))S( Vl - V(t n+1 ; (xi, Vj ),t n ))\ 

k,l 

< iXit^-ix^v,),^) - X h (t n+1 ;( Xi , Vj ),t n )\, 

and 

k,l 

< |y(t" +1 ;( 2;i ,t; J ),t n )-T4(t" +1 ;(x i ,^),t n )|. 

So that we get: 

\\Rhfh +1 - Rhfh +1 \\Ll(Q) < AxA vJ2\ui,j{fh) n \(\X(t n+1 ; (xi,Vj),t n ) - X h (t n+1 ;(xi,Vj),t n )\ 

+ \V(t n+1 ;(x t ,v j ),t n )-V h (t n+1 ;(x l ,v j ),t n ))\), 
^max^^ 1 ;^),^^^ 

thanks to particle conservation (Yli j ^ijifh) = Y2i j w «j(/o)) and positivity of Uij(f^), where 
(X,V) are the exact characteristics, solution of the differential system (|2.7p . and (Xh,Vh) the 
approximate characteristics defined in (13. lip . 

To move on, we need another lemma which enables to control the difference between exact and 
computed characteristics. It clearly depends on the algorithm we use. Let us first give the lemma 
for the Verlet algorithm. 

Lemma 7 : Verlet 

If E G W 2 '°°([0, t] x K), and with (X,V) calculated exactly with the differential system <\2.7h . 
and (Xh, Vh) computed with and a Verlet algorithm: 

\x{t n+i -{ XuVj ),n x^r+ i ;(x i) ^o ) t n )i + in* n+1 ;(^,^),* n )-^(t n+ ^(^,^),ni 

< C^ + kt\\{E-E h ){t n+1 2)\\ L ™. 

Proof: The strategy follows the work of M. Bostan and N. Crouseilles ([2]). 

Let us recall how (x,v) n+1 is computed from (x,v) n with the Verlet algorithm 

x n + \ = x n + 

v n+l =v n + AtE h (t n+ ^,X n + — V n ), 
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x n+l = x n+± + ^ v n+\ 

Then we define X h (t n+1 , (x n ,v n ),t n ) = x n+l and V h (t n+1 , (x n ,v n ),t n ) = v n+l . 
Let us begin with the characteristics in v: 

(V h -V)(t n+1 ;(x n ,v n ),t n ) = v n + AtE h (t n+1 2,x n + —v n )-v n - E(s, X(s; (x n ,v n ),t n ) ds, 

2 .l+n 



(E(s, X(s; (x n ,v n ),t n ) - E(t n+1 * , X{t n+ h- ( x n ,v n ),t n ))) ds, 
At(E(t n+ ^,X(t n ^; (x n ,v n ),t n )) - E h (t n+ ^,x n + — v n )). (4.29) 



Let us take care of the integral term in (|4.29p . using a Taylor expansion around s = t n+ 2 of 
s i-> E{s,X{s;{x n ,v n ),t n ). 

E(s, X(s; (x n ,v n ),t n ) = E(t n+ ^,X(t n+ ^;(x n ,v n ),t n )) + (s - t n+ ^)E' (t n+ ^ , X(t n+ ^ ; (x n ,v n ),t n )) 

+ / As-u)E"{u,X{u;{x n ,v n ),t n ))du. (4.30) 

Jt n+ ? 

Let us precise that E'(s,X(s)) = £e(s,X(s)), E"(s,X(s)) = £E'(s,X(s)). Then using (PUD in 
(|429jh we get 



n+l 

{E(s,X(s;(x n ,v n ),t n )) - E(t n+ ^,X{t n+ ^,(x n ,v n ),t n ))ds, 

= E'{t n ^ , X(t^\ ; (x» v n ), t n )) [(i^i) 2 ]*: + 1 , 

+ / / (s-u)E"(u,X{u;(x,v),t n )duds. 
Jt n Jt n+ ? 

There are here two terms to control. The first one is zero, and for the second one, we have, 

t n+1 ps rt n+1 t-s 

As -u)E"{u,X{u;{x n ,v n ),t n )duds\ < \\E"\\ L ~> / (s-u)duds, 

t n+l 



< \\e"\\l~ J tn [-(^) 2 i; + ,&, 

< f (s-t n ^) 2 ds, 

< - w \\E"\\ L o C <C\\E"\\ Lac At 3 . (4.31) 



Now let us deal with the second term of (|4.29p . Since E is bounded: 



\x n + ^v n -X(t n+L 2-(x n ,v n ),t n )\ = | / v n -V{s;(x n ,v n ),t n )ds\, 

< / (s-t n )V(u;(x n ,v n ),t n )ds u£[t n ,s], 



< C(\\E\\ L ~) [ {s - t n ) ds < C At 2 . (4.32) 
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and thus, since E' is bounded, using the zero mean theorem: 

\E(t n+ ^,X(t n+ h-( x n ,v n ),t n )) - E(t n+ h, x n+ ^)\ < C(\\E'\\ L o )\X(t n+ h;(x n ,v n ),t n ) - x n+ ?\, 

< C{\\E'\\ L ~)At 2 . (4.33) 

Finally the second term of (I4.29P can be controlled by 

\E(t n+1 i,X(t n+1 2;(x n ,v n ),t n )) - E h {t n+1 *,x n+l 2)\ < \\(E - E h )(t^)\\oo 

+ \E(t n+ ^,X(t n+1 2-(x n ,v n ),t n )) - E(t n+ Kx n+ ^)\. (4.34) 

So, using (|PT|) . (031) and (Olj) . we get 

104 - V)(t n+1 ; (x, v), t n )\ < CAt 3 + At\\(E - E h )(t n+L 2)\\ Lao(n) . (4.35) 
Let us now deal with the characteristics in X: 



At „ At 



t n+l 



(X h -X)(t n+1 ;(x n ,v n ),t n ) = v n + —v n+l - V(s;(x n ,v n ),t n )ds, 

Z I Jpi 

= -At(V(t n+ ^;{x n ,v n ),t n ) - ~v n - ^v n+1 ), 

t n+l 

(V(s; (x n , v n ),t n ) - V(t n+1 z; (x n ,v n ),t n ))) ds, (4.36) 

i 

so once again we have to control two terms. 

For the first one, thanks to Taylor's inequality, like for X, it comes: 

(V(s;(x n ,v n ),t n ) -V{t n+1 2;{x n ,v n ),t n ))ds\ < C(||£|| Loo )At 3 . (4.37) 

Let us precise that this is nothing else than the error in the mid-point rule for numerical integration. 
Now, the second term in (|4.36|i : 

t n+ i 



V(t n+ 2-( x n ,v n ),t n )) = v n + E(s,X(s;(x n ,v n ),t n ))ds, 

= v n + ^E(t n+1 2, X{t n+1 2- (x n ,v n ),t n )), 



and 



rt «+4 

+ / (E(s,X(s; (x n ,v n ),t n )) -E(t n+ 2,X(t n+ 2;{x n ,v n ),t n )))ds, 



(E(s, X(s; (x n ,v n ),t n )) - E{t n+l 2,X{t n+l 2 ; (x n , v n ), t n ))) ds\ < C(||£'|| LC o)At 2 



with the error formula for the rectangle rule. On the other hand 
V(t n+ *;{x n ,v n ),t n )-±v n -\v n+1 = l -v n - l -v n+1 



+ ^E(t n+1 * , X(t n+L 2- (x n ,v n ),t n )) + 0(At 2 ). 
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Since 

i( u «_ t ,«+l) = -^ h (t n+ 3,X n+ 3), 

we have, proceeding as for (|4.33j) - (|4.34p 

\V(t n+ ^; (x n ,v n ),t n ) - ~v n - \v n+1 \ = \^{E(t n+1 2,X(t n+ ^;(x n ,v n ),t n )) - E h {t n+ ^,x n+ ^) 



+0(At 



2\ 



At 



< —(\\(E-E h )(t n ^)\\L^(n) + CAt'). (4.38) 

To conclude, using (|4.37p . (|4.38p . we have: 

\(X h - X)(t n+1 ; (x,v),t n )\ < CAt 3 + CAt 2 \\(E - E h )(t n+1 2)\\ LOO{n) . (4.39) 

Finally, using (|4.35p and (|4.39p . we get the estimation of Lemma 7, and using (|4.28p . this also 
implies 



\\R h fk n+1 - RhfZ +1 \\ L u a) < CAt 3 + At\\(E-E h )(t n +^)\\ L 



oo . 



Lemma 8 : CK3 

If E E W 4 '°°([0, t] x M), and with (X,V) calculated exactly with the differential system (|2.7j) . 
and (Xh,Vh) computed with Eh,Ph,Jh and a CK3 algorithm: 

\x{t n+l -(x h v 3 ),t n ) x fe (r+ 1 ;(x i ,^o,ni + in* n+1 ;(^,^),n-^(* n+1 ;(^,^),i n )l 

< CAt 4 + C{At\\{E n - ^)||,- (n + At 2 \\(r - fflWl-M 
+At 3 ||(^-<tf)||,« (n) ). 

Proof: This proof just relies on Taylor expansions and computations already made: 

At 2 At 3 

X(t n+1 ; (xi, Vj ),t n ) = Xi + Atvj + — E n ( Xi ) + —^( Xi)Vj ) + 0(At 4 ), 

At 2 At 3 
X h (t n+1 ; (xi, Vj ),t n ) = Xi + Atvj + —E%( Xi ) + — ffifruVj), 

and 

At 2 At 3 

V(t n+1 ;( Xi ,v 3 ),t n ) = vj + AtE n (xi) + —<p n ( Xi , Vj ) + —<p n (xi,Vj) + 0(At 4 ), 

At 2 At 3 
V h {t n+1 ; ( Xi , Vj ),t n ) = Vj + AtEJt(xi) + —ftixiw) + —rtfavj), 

and the lemma follows by simple subtraction. 

In both cases we need to control the difference between the exact and approximate fields. Let 
us begin with Verlet algorithm. 

Lemma 9 : Verlet 

If EG. W 2 '°°([0,t] x R) , it comes 

\\(E- E h )(t n+1 2)\\ L °°(n) < C{h 2 + At 2 + Ath 2 + e n ) 
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Proof: First 

E(t n+1 2,x)= [ K(x,y)(l f(t n+ ^,y,v)dv-l)dy, 



o 

rL 

E h (t n +i,x)= I K[x,y){\ f h (t n+± 2,y,v)dv-l)dy. 

Hence 



o 



E(t n+ i,x)-E h (t n+ i,x) = t ' K(x,y)U {f{t n+ i,y,v)-f h (t n+ i,y,v))dv)dy, 



o 



L 

n+ 











K(x, y){ / (f(t n+ 3 ,y,v)- R h f(t n+ 2 ,y,v)) dv) dy, 
Jr 

K(x, y)( / (R h f(t n+ 2 , y, v)- f(t n+ z , y, v)) dv) dy, 
Jr 

K(x, y){ [ (f(t n+l 2 , y , v )- f h (f + \ , y , v )) dv) dy, 
Jr 

K(x,y)([ (f h (t n+1 2,y,v)-f h (t n+1 2,y,v))dv)dy, (4.40) 



o 



where 



f(t n+ Ky,v) = J2^Af n )S(y - X(f+h (x k , Vl ),t n )S{v - V(t n+1 *; (x k , Vl ),f 



k.i 



and 



, y, v) = J2 »>M)S(V ~ X(t n+l ; (a*, vi), t n )S(v - V(t n+ h ; (x k , Vl ),t n ). 
k,l 



In order to lighten notations, X(t n+ z ; (x k , vf), t n ) and V(t n+ 2 • (x k , vi), t n ) will be denoted X k ~^ 2 

and V k ■ 2 . We have four terms to control. 

The first one is controlled using property (ii) of consistency and accuracy: 

L 



Now, the second term of (|4.40|) . 

L K(x, y){ [ (R h f(t n+1 i ,y,v)- f(t n+L 2 , y , v )) dv) dy, 
Jr 



K{x,y){ {f(t n+ 2,y,v)-R h f{t n+ i,y,v))dv)dy\<C\\K\\ 00 h\ (4.41) 



o 



K(x, y)( / (E(u%[*S(y - x k )S(v - «,) - ^(y - X n k J*)S(v - V™ +5 ))) dv) dy. 



k,l 

We have, using Taylor expansion and Vlasov equation: 



= f n+ hx k ,v l ) = f n (x k ,v l ) + ^^(x k ,v l )+0(At 2 ), 

= P(xk,vi) - ^(v l ^(x k ,v l ) + E n (x k )^(x k ,v l )) + 0(At 2 ) 
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Moreover, since S is piecewise polynomial of degree one and continuous, we have almost everywhere 
(which is enough as we are going to integrate these expressions) 

S(y - x£*) = S(y -x k - ^ + 0(At 2 )) = S(y - x k ) - ^■v l S'(y - x k ) + 0(At 2 ). 
Using, these expansions, we get: 

i-L 



[ K(x, y)( [ ($>Jt*Sfo ~ Xk ^ v ~ ^ ~ u h S ^ ~ X kf 2 )S(v - V^ 1 ))) dv) dy, 

At f L f x df n 

= ir[ / K (x,y)( (S^viSiv - vi)(f n (x kl vi)S'(y - x k ) - — — (x k ,vi)S(y - x k )), 

^ M ki 

+ E n {x k )S{y - x k ){f n (x k ,v l )S'{v - v t ) - ^(x k , v t )S(v - v t ))) dv) dy) + 0(At 2 ). (4.42) 
k,l ° X 

There are two terms in (j4.42p . They will be dealt with similarly using mid-point quadrature, 
which is of second order. For the first term, it writes 

K (x,y)( / (y]v t S(v - vi)(f n (x k ,vi)S'(y - x k ) - -^—(x k ,vi)S(y - x k ))dvdy, 

= Av Y] vi / K(x, y)(f n (x k ,vi)S'(y - x k ) - ——(x k ,vi)S(y - x k )) dy, 
k,l Jo dx 

= Av^2(viK(x,x i+ i)(f n (x k ,vi)S'(x i+ i - x k ) - -Q—(x k ,vi)S(x i+ i -x k )). 

i,k,l 

Here, we have to use the properties of linear splines. S'(x i+ i — x k ) and S(x i+ i — x k ) are non 
zero only if k = i or k = i + 1. Then, we have: S'(x i+ i — Xi) = S'(x i+ i — Xi + \) = ^ and 
S(x i+ i — Xi) = S(x i+ i — Xi + \) = |. Using that, we get 

Av S ^{viK(x,x i+ i)(f n (x kl vi)S'{x i+ i - x k ) - -^—{x k ,vi)S(x i+ i - x k )), 

2 2 CJX 2 

i,k,l 

= Av^ Vl K{x,x i+ i) I — l—{x i ,v l ) + —(x i+1 ,v l ) 

i,l V V 

Using Taylor expansions with respect to x, we easily get: 

/(x, +1 ,^-/(x„„) = ^ (Xi+ ,, [)+0(Ax2) . 

and 

1 f df n df n \ df n 

2 \~dx~( XuVl) + ~dx~^ Xi+l ' Vl) ) = ~dx~ i ~ Xi+1 2 }Vl) + °( Ax2 )' 

since / G Wc^ rx (R x x R v ). And to conclude for the first term of (|Q2]> : 

•L r Qfn 
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K(x,y)(f (Y i viS(v-v l )(f n (xk,vi)S'(y-x k ) 
7k u 



(x k ,vi)S(y - x k ))dvdy\ < CAx . 
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For the second term of (I4.42p . with a mid-point quadrature for the integral with respect to v, and 
the same properties of splines: 

K(x, y)( / (V E n {x k )S{y - x k )(f n (x k , Vl )S'(v - Vl ) - -£-(x k) vi)S(v - «,))) dv) dy, 
k,l dx 

= Av^E n (x k ) / K(x,y)S(y - x k ) dy ^{f n {x k , Vi)S'(v j+ i - vi) - -Q—(x k ,vi)S(v j+ i - vi)), 
k.i Jo j X 

< AxAv\\K\\ L oo y \E n (x k )\\ J - — J - -(— — (x k ,v j+1 ) + -5— (x k ,Vj)) \. 

' 1 At; 2 ov ov 

Using again Taylors expansions, with respect to v, we get: 

L K(x, y)( / (V E n (x k )S(y - x k )(f n (x k , Vl )S'(v - «,) - ^f-(x k , Vl )S(v - v t ))) dv) dy, \ 
M 

< AxAv\\K\\ Laa Y,\E n (x k )\Av 2 , 

k,j 

<C\\E n \\ Loo{[QM) Av 2 <C'Av 2 , 
since E G iy 2,oo ([0, t] x R). To conclude, the second term of ([4.40p can be bounded like that: 

L K(x,y)([ (R h f(t n+ Ky,v)-f(t n+ Ky,v))dv)dy\<C(Ath 2 + At 2 ). (4.43) 



o 



o 



For the third term of (|4.40[) : 



K(x, y)( I (f(t n+ 2 ,y,v)- f h (t n+ 2, y, v)) dv) dy, 
-L 



= | / K(x,y)(! Y,Mf n )-^Afk))S(y-X:^)S(v-V k n ^)))dv)dy, 
Jo Jr kil 

<AvY" MD ~ "k,lUh)\ [ L \K(x, y)\S(y - xfi*)dy, 
k,i Jo 

< WKWooAxAvY.iri^vi) - fj:(x k , Vl )\ <Ce n . (4.44) 
k,l 

Eventually, the last term of (|4.4Qj) : 

L K(x, y)( [ (f h (t n+l 2 ,y,v)- f h (t n+l 2 , y, v)) dv) dy\, 







= 1 / K{x,y){f Y,^M(s(y-x n k ^hs(v-v^h-s(y-x^ 

< Av^M) [ L K(x,y)(S(y - X^) - S(y - dy, 

k~l Jo 

< \\K\\^AxAv^M)\KP - X tkh 

k,l 

< CAt 2 , (4.45) 
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using successively the positivity of fh, mass conservation, fh is 1-lipschitzian and a result in Lemma 

7: (|4.32p . where obviously, i^h-ikiy^h-Ckl)^ are ^ e a PP ox i ma te characteristic curves at time t n+ 2 
beginning at time t n at (xk,vi). 

To conclude, using (itlTi (14431) . (Oil) and (EU5jl we get: 

IK-E- ^)(t n+ 5)||L-(Q) < C(/* 2 + At 2 + e n + Ath 2 ), 
which is what was expected. 

Lemma 10 CK3 With the same hypothesis as in Lemma 8, we have : 



A+ 3 A-h 2 h 2 h 2 

\\R h f h n+1 -R h f h n+1 \\ L}{n) < C(e n (At + ^L + ^-)+Ath 2 + At 2 ^- + At^ + At*). 

hK ' Ax z Ax Ax Ax z 



Proof: Here, we need to evaluate the difference between the £°°(Q) norms of the exact and ap- 
proximate values of (f>, <p and tp, so to say the one between: 

• E n and Eft, p n and pJJ, J n and J™, 

• their first spatial derivative and the one of IV; and IV^ . 
Let us start with 



(E n -El){ Xl ) = [ K{x u y){f 

JO JR 

< C(e n + h 2 ). 



(f(t n ,y,v)-f h (t n ,y,v)dv)dy, 



simply using a quadrature with the mesh points, which will also be of second order thanks to 
periodicity, and the fact that K is bounded. So that: 

\\E n -E n h \\ laa{%L]) <C{e n + h 2 ), (4.46) 
\lP n - P n h\\LU[o,L]) = Ax Y,\ I (f U ~ fh)Mdv\, 

i JR 

< C(e n + h 2 ). 

So that using the equivalence of discrete norms, carefully noticing that 

1 ,, 



we get: 



\P ~ Ph\\l°°([0,L]) ^ ~ Ph\\Ll([0,L])' 

c 



and 



n , z,2\ 



< C(e n + h z ), 
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using the same arguments and the fact that f is compactly supported, so that: 

1 1 in jn 1 1 ^ 1 1 jn jn 1 1 

\W - J h\\l™([0,L]) - £^W J ~ J h\\Ll([0,L})i 

Let us precise that the same bound is obviously also valid for J. So that we get, still using that f2 
is bounded: 

ll^-€ll^([o,L])<^(e n + / i 2 ). 

For the first spatial derivative of these three terms, we can use the same strategy of finite 
difference. Let us do it with E. 

d(E - E h ) _ (gVi - ar+i,J - (gLi - E?_ lih ) 

dx 2Ax +U{IXx ), 

so that using (|4.46p we get: 



d(E - E h ) C 2 

1 dx {t ~ Ax- ( 



-n\\i~(io, L] )<^ n +n 

For p n , J n and for IV; just bounding v, v 2 in its integral definition, the same strategy leads to 

ll^-^ll^([o,L])<^2(e n + / i 2 ). 
Plugging these estimates into Lemma 8 and then Lemma 6 completes the proof. 

4.3 End of the proof 

For the sake of simplicity, and since we are interested in At, Ax tend to 0, we will assume At < 1 
We can now apply Gronwall inequality since 

Verlet For Verlet, using Lemmas 4, 5, 7, 9, we get, 

e n+1 < C{h 2 + Ai 3 + (At 2 + At){h 2 + h 2 AtAt 2 + e n ) + e n , 

< (1 + CAt)e n + C(h 2 + At 3 + (At 2 + At)(h 2 + At 2 )), 

So that 

h 2 

e n < exp(C / T)e° + C(h 2 + At 2 + — ), 
which is what was expected. Taking At = Ch a , we find the best global order with a = | being |. 

CK2 For CK2, using Lemmas 4, 5, 8, 10, we get 

At 2 

e n+1 < C(h 2 + At 3 + (— + At)(h 2 + e n )) + e n , 

Ax 

At 2 At 2 

< (1 + C(At + =-)e n + C(h 2 + At 3 + At h 2 + ^-h 2 ), 

Ax Ax 

Here, assuming At < Ax, we have : 

e n < exp(C'T)e° + C(h 2 + At 2 + ^), 

If you want to look for the best global order here, you find the same result as in Verlet, nevertheless, 
this cannot fit with the other assumption At < Ax. Therefore, the only way is to take At = Ax, 
and the global order is 1. 
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CK3 For CK2>, using the same lemmas as for CK2: 

e n+i < C (h 2 + At 4 + e n (At + ^ + ^) + Ath 2 + At 2 ^ + At^) + e n , 

Ax Ax A Ax Ar 

< e"(l + C(At + ^ + + °( h2 + Ath2 + + + Ai ')- 

Ax Ax z Ax Ar 

Assuming again At < Ax, we have 

h 2 

e n < exp(C'T)e° + C(At 3 + h 2 + ^). 

The same remark as with CK2 is still valid. We can see that we are limited because of the terms 
Si anc ^ sf ■ ^ n or< ier to be able to reach higher orders, we would have to use splines of superior 
degrees m > 1 to get terms in like in the other proofs of convergence, for example: ([U[5]) 

5 Numerical results 

In order to validate our new schemes, we have tested them on two standard test cases of plasma 
physics, the two stream instability and the bump on tail instability. We also compared them to the 
classical and knowledgeably robust Verlet scheme. Notice that because of the diffusitivity of linear 
splines, we have used cubic splines for the distribution function. 
For the two stream instability, the initial condition is given by 

1 2 

fo(x,v) = e~" > 2 v 2 [l — acos(kx)], 

V 2ir 



with k = 0.2 and a = 0.05. The computational domain is [0,2-7r/A;] x [—9,9] which is sampled 
by N x = N v = 128 points. We used a time step At = 0.1 in the results of the left hand side of 
Figure [T] and of At = 0.3 on the right-hand side of the figure. We display the L 2 norm which 
reveals the dissipation of the scheme, the total momentum and the total energy. All of those are 
conserved in the continuous Vlasov-Poisson system. We do not display the number of particles 
which is conserved with an even better accuracy than the momentum. The momentum is exactly 
conserved by the scheme and up to about 10~ 13 in the simulation. This is due to roundoff errors 
and the truncation of the velocity space. The L 2 norm cannot be exactly conserved by any scheme 
using a phase space grid as soon as the grid does not resolve anymore the filaments. The Verlet 
scheme is our reference scheme here, and we observe that the results obtained with the CK schemes 
are very close, especially for the smallest time step. Moreover conservation properties are better 
for the third order CK3 than for the second order CK2. 

For the bump-on-tail instability test case the initial condition writes 

fo{x,v) = f(v)[l +acos(kx)], 

with 



f(v) = n p exp(— v 2 /2) + n^exp 



I 1 2 

\v — u\ 



2v 2 

on the interval [0, 20-7r], with periodic conditions in space. The initial condition fo is a Maxwellian 
distribution function which has a bump on the Maxwell distribution tail; the parameters of this 
bump are the following 

9 2 

1T-V = ; --7777 , Tlh = : 7T777i U = 4.5, Vf = 0.5, 

p 10(2vr) 1 /2' 10(2tt) 1 /2' 
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Figure 2: Bump on tail 



whereas the numerical parameters are N x = 128, N v = 128, v m ax = 9, At = 0.2. The results are 
displayed in Figure [2j Here the momentum does not vanish, so that the results are not polluted 
by roundoff errors and the momentum is exactly conserved. The conclusion for the L 2 norm and 
the total energy is the same as in the Two Stream Instability test case. The potential or electric 
energy is a classical diagnostic for the bump on tail instability. The oscillations go on for a long 
time with all three time schemes, even though there is a slight energy increase for the CK2 scheme. 



6 Conclusion 

In this paper, the proof of a L 1 convergence has been reached for linear spline interpolation. The 
originality, except from the choice of the L 1 norm is that the convergence has been reached for a 
non split method. In this paper, the computation of the characteristics has been made with the 
Verlet algorithm, or with a CK procedure, but the proof can be adapted to other algorithms such 
as Runge Kutta of any order. There remains for the moment some problems using splines of higher 
orders, especially concerning stability. This prevents us from reaching real high order algorithms. 
Numerical experiments that can be seen in [7j, and confirmed here seem to prove that the method 
is also stable and convergent for cubic splines. Nevertheless, there remains a problem to preserve 
the I 1 norm of the coefficients OJij, since some of them can become non positive in the solving of 
the linear system with splines of degree higher than 2. Another way of tackling the problem will 
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probably be needed. 
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